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Abstract 

The focus of this work is on rigorous mathematical analysis of the topological derivative 
based detection algorithms for the localization of an elastic inclusion of vanishing charac- 
teristic size. A filtered quadratic misfit is considered and the performance of the topological 
derivative imaging functional resulting therefrom is analyzed. Our analysis reveals that the 
imaging functional may not attain its maximum at the location of the inclusion. Moreover, 
the resolution of the image is below the diffraction limit. Both phenomena are due to the 
coupling of pressure and shear waves propagating with different wave speeds and polariza- 
tion directions. A novel imaging functional based on the weighted Helmholtz decomposition 
of the topological derivative is, therefore, introduced. It is thereby substantiated that the 
maximum of the imaging functional is attained at the location of the inclusion and the 
resolution is enhanced and it proves to be the diffraction limit. Finally, we investigate the 
stability of the proposed imaging functionals with respect to measurement and medium 
noises. 
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1 Introduction 

We consider the inverse problem of identifying the location of a small elastic inclusion in a ho- 
mogeneous isotropic background medium from boundary measurements. The main motivations 
of this work are Non-Destructive Testing (NDT) of elastic structures for material impurities 
[P5] . exploration geophysics [T], and medical diagnosis, in particular, for detection of potential 
tumors of diminishing size |25j . 
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The long standing problem of anomaly detection has been addressed using a variety of 
techniques including small volume expansion methods [8j [9] , MUSIC type algorithms [4] and 
time-reversal techniques [5J [5] . The focus of the present study is on the topological derivative 
based anomaly detection algorithms for elasticity. As shown in [5], in anti-plane elasticity, the 
topological derivative based imaging functional performs well and is robust with respect to 
noise and sparse or limited view measurements. The objective of this work is to extend this 
concept to the general case of linear isotropic elasticity. The analysis is much more delicate in 
the general case than in the scalar case because of the coupling between the shear and pressure 
waves. 

The concept of topological derivative (TD), initially proposed for shape optimization in 
[T51 [2H U2]> has been recently applied to the imaging of small anomalies, see for instance, 
[131 [Ml IT7l [TBI [19l [20l [23] and references therein. However, its use in the context of imaging 
has been heuristic and lacks mathematical justifications, notwithstanding its usefulness. 

In a prior work [5] , acoustic anomaly detection algorithms based on the concept of TD arc 
analyzed and their performance is compared with different detection techniques. Moreover, a 
stability and resolution analysis is carried out in the presence of measurement and medium 
noises. 

The aim of this work is to analyze the ability of the TD based sensitivity framework for 
detecting elastic inclusions of vanishing characteristic size. Precisely, our goal is threefold: 
(i) to perform a rigorous mathematical analysis of the TD based imaging; (2) to design a 
modified imaging framework based on the analysis. In the case of a density contrast, the 
modified framework yields a topological derivative based imaging functional, i.e., deriving from 
the topological derivative of a discrepancy functional. However, in the case where the Lame 
coefficients of the small inclusion arc different from those of the background medium, the 
modified functional is rather of a Kirchhoff type. It is based on the correlations between, 
separately, the shear and compressional parts of the backpropagation of the data and those of 
the background solution. It can not be derived as the topological derivative of a discrepancy 
functional; and (3) to investigate the stability of the proposed imaging functionals with respect 
to measurement and medium noises. 

In order to put this work in a proper context, we emphasize some of its significant achieve- 
ments. A trial inclusion is created in the background medium at a given search location. Then, 
a discrepancy functional is considered (c.f. Section [3]), which is the elastic counterpart of the fil- 
tered quadratic misfit proposed in [5] . The search points that minimize the discrepancy between 
measured data and the fitted data are then sought for. In order to find its minima, the misfit is 
expanded using the asymptotic expansions due to the perturbation of the displacement field in 
the presence of an inclusion versus its characteristic size. The first order term in the expansion 
is then referred to as TD of the misfit (c.f. Section l3~T|) which synthesizes its sensitivity relative 
to the insertion of an inclusion at a given search location. We show that its maximum, which 
corresponds to the point at which the insertion of the inclusion maximally decreases the misfit, 
may not be at the location of the true inclusion (c.f. Section I3~l2"j) . Further, it is revealed that 
its resolution is low due to the coupling of pressure and shear wave modes having different wave 
speeds and polarization directions. Nevertheless, the coupling terms responsible for this degen- 
eracy can be canceled out using a modified imaging framework. A weighed imaging functional 
is defined using the concept of a weighted Helmholtz decomposition, initially proposed in [3] 
for time reversal imaging of extended elastic sources. It is proved that the modified detection 
algorithm provides a resolution limit of the order of half a wavelength, indeed, as the new func- 
tional behaves as the square of the imaginary part of a pressure or shear Green function (c.f. 
Section ET2"1) . For simplicity, we restrict ourselves to the study of two particular situations when 
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we have only a density contrast or an elasticity contrast. In order to cater to various applica- 
tions, we provide explicit results for the canonical cases of circular and spherical inclusions. It 
is also important to note that the formulae of the TD based functionals arc explicit in terms of 
the incident wave and the free space fundamental solution instead of the Green function in the 
bounded domain with imposed boundary conditions. This is in contrast with the prior results, 
see for instance, [TS]. Albeit a Neumann boundary condition is imposed on the displacement 
field, the results of this paper extend to the problem with Dirichlct boundary conditions. A 
stability analysis of the TD based imaging functionals was also missing in the literature. In 
this paper we carry out a detailed stability analysis of the proposed imaging functionals with 
respect to both measurement and medium noises. 

The rest of this paper is organized as follows: In Section [21 we introduce some notation 
and present the asymptotic expansions due to the perturbation of the displacement field in 
the presence of small inclusions. Section [3] is devoted to the study of TD imaging functional 
resulting from the expansion of the filtered quadratic misfit with respect to the size of the 
inclusion. As discussed in Section 13.21 the resolution in TD imaging framework is not optimal. 
Therefore, a modified imaging framework is established in Section 01 The sensitivity analysis 
of the modified framework is presented in Section 14.21 Sections [5] and [5] arc devoted to the 
stability analysis with respect to measurement and medium noises, respectively. The paper is 
concluded in Section [7] 

2 Mathematical formulation 

This section is devoted to preliminaries, notation and assumptions used in rest of this paper. 
We also recall a few fundamental results related to small volume asymptotic expansions of the 
displacement field due to the presence of a penetrable inclusion with respect to the size of the 
inclusion, which will be essential in the sequel. 

2.1 Preliminaries and Notations 

Consider a homogeneous isotropic elastic material occupying a bounded domain il C M. d , for 
d = 2 or 3, with connected Lipschitz boundary dft. Let the Lame (compressional and shear) 
parameters of f2 be Ao and po (respectively) in the absence of any inclusion and pa > be the 
(constant) volume density of the background. Let D C be an elastic inclusion with Lame 
parameters Ai, /ii and density p\ > 0. Suppose that D is given by 



where B is a bounded Lipschitz domain in R containing the origin and z a represents the 
location of the inclusion D. The small parameter 6 represents the characteristic size of the 
diameter of D. Moreover, we assume that D is separated apart from the boundary dfl, i.e., 
there exists a constant Co > such that 



D := SB + z 



a 



(2.1) 



inf dist(x, 90) > Co, 



(2.2) 



where dist denotes the distance. Further, it is assumed that 



d\ m + 2^ rn > 0, n m > 0, m e {0, 1}, (Ao - Ai)(^ - Mi) > 0. 



(2.3) 
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Consider the following transmission problem with the Neumann boundary condition: 



u _ = u| 

du 

dv 



p UJ 2 U 

piw 2 u 



I — 

du 

dv 



du 

dv 



in 0\D, 
in D, 
on <9£>, 

on 3D, 
on 



(2.4) 



where u > is the angular frequency of the mechanical oscillations, the linear elasticity system 

d 

£\ no an d the co-normal derivative — , associated with parameters (Ao, Mo) are defined by 

dv 

(2.5) 



£a ,mo I w ] : = MoAw + (A + Mo)VV • w 



and 



dw 

dv 



A (V • w)n + ^o(Vw J + (Vw J Y )n, 



(2.6) 



respectively. Here superscript T indicates the transpose of a matrix, n represents the outward 
unit normal to dD, and J= is the co-normal derivative associated with (Ai,Mi)- To insure 
well-posedness, we assume that pqoj 2 is different from the Neumann eigenvalues of the operator 
mo m (^ 2 (^)) ■ Using the theory of collectively compact operators (see, for instance, 
Appendix A.3]), one can show that for small <5 the transmission problem (|2.4p has a unique 

solution for any g e (L 2 (d£l)) d . 

Throughout this work, for a domain X, notations |_ and |+ indicate respectively the limits 
from inside and from outside X to its boundary dX, Sij represents the Kroneckcr's symbol and 



a,(3 e {P,S}, i,j,k,l,i',j',k',l',p,q£ {1, • • • ,d}, 
where P and S stand for pressure and shear parts, respectively. 



e{o,i}, 



Statement of the Problem: 

The problem under consideration is the following: 

Given the displacement field u, the solution of the Neumann problem (|2.4|) at the boundary 
dQ, identify the location z Q of the inclusion D using a TD based sensitivity framework. 



2.2 Asymptotic analysis and fundamental results 

Consider the fundamental solution Pj^(x, y) := T^(x — y) of the homogeneous time-harmonic 
elastic wave equation in K d with parameters (A m , p m , p m ), i.e., the solution to 

(£\ m ,» m + PmW 2 )^(x - y) = d y (x)I 2 , Vx £l d ,x/ y, (2.7) 
subject to the Kupradze's outgoing radiation conditions |22j . where 8 y is the Dirac mass at y 
and I 2 is the d x d identity matrix. Let c s ~ ^ anA - / Ao+2 ^° 



and c P 
po V P" 

and the pressure wave speeds respectively. Then Tq is given by [T] 



be the background shear 



I 2 Gg(x) 



1 

p UJ 2 



[G£(x) - Gg(x)] 



x € 



d = 2,3, 



(2. 
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where the tensor D x is defined by 

D x = V x ® V x = (%)f J=1 , 
and the function is the fundamental solution to the Hchnholtz operator, i.e., 

(A + 4)G£(x) =<5 (x) xel^.x/O, 
subject to the Sommerfeld's outgoing radiation condition 

■' (x) =o{R 1 - d / 2 ), x£dB{0,R), 



with -6(0, R) being the sphere of radius R and center the origin. Here d. 

is the wave-number, and ^ represents the normal derivative. 
The function GJ^ is given by 



-IffW^lxl), d = 2, 



d = 3, 



(2.9) 



471-1x1 ' 



where Hn' is the order n Hankel function of first kind. 

Note that Tq can be decomposed into shear and pressure components i.e. 

r-(x) =r£ s (x) + r£ P (x), VxeM d , x^o, 

where 



1 



1 



r£ P (x) = l T D x G^(x) and r£ s (x) = -^-(«|l 2 + D x )Gg(x). 

;e that V • T^ s = and V x T% p = 0. 

Let us define the single layer potential S£ associated with (£a .^ + Po^ 2 ) by 



(2.10) 
(2.11) 



SS[*](x) := f r-(x-y)*(y)d<7(y), 



x e 



and the boundary integral operator Kf, by 



^[*](x) 



p.v. 



/ — rS(x-y)#(y)dtr(y), a.e. x e Sf2 



(2-12) 



(2.13) 



for any function $ € (L 2 (<9f2)) d , where p.v. stands for Cauchy principle value. 



Let (1Cq)* be the adjoint operator of on (L 2 (dQ)) , i.e., 

(iC?J)*[*](x) = p.v. f /-r-(x - y)*(y)My), a.e. x e 



dfl. 



It is well known, see for instance [21 Section 3.4.3], that the single layer potential, Sq, enjoys 
the following jump conditions: 



g(gB[*l) 

3^ 



(x)= ±- J [S](x), a.e. xean. 



(2.14) 
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Let N w (x, y), for all y € SI, be the Neumann solution associated with (Ao, po, Po) m *- e -, 

(£ Ao ,mo + /W)N w (x,y) = -5 y {x)I 2 , x G Q, x ^ y, 
<9N^ 



9i> 



(x,y)=0 



x e on. 



(2.15) 



Then, by slightly modifying the proof for the case of zero frequency in [5], one can show that 
the following result holds. 



Lemma 2.1. For all x 6 <9£1 and y £ Q, we have 

1 



For i,j e {1, 



7 + /C£ [iSr(.,y)](x)=r£(x-y) 



(2.16) 



, d}, let v,j be the solution to 

Cxa.fioVij = 
£Ai,/*iVij = 



9v, : 



y 1 + 



in lr\5, 
in 5, 
on 

on 95, 



(2.17) 



v,:j(x) — Xiej — O (|x| 1 d ) as |x| — » oo, 

where (ei, • • • , e<j) denotes the standard basis for M. d . Then the elastic moment tensor (EMT) 
M := ( m ijpq)jj pq -i associated with domain B and the Lame parameters (Ao, /Uoj ^l, Ml) is 
defined by 

~d(x p e q ) d{x p e q ) 



Hjpq 



OB 



Op 



dv 



Vij da, 



(2.18) 



can be expressed as 



see [81 111). In particular, for a circular or a spherical inclusion, 

M = aI 4 + 6I 2 ®l2, (2.19) 

or equivalently as 

rriijki = - (SikSji + 5ii5 jk ) + bSySki, 

for some constants a and b depending only on Ao,Ai,^o,/ii and the space dimension d [2 
Section 7.3.2]. Here I4 is the identity 4-tensor. Note that for any d x d symmetric matrix A, 
14(A) = A. Furthermore, throughout this paper we make the assumption that fix > and 
Ai > Ao in order to insure that the constants a and b are positive. 
It is known that EMT, M, has the following symmetry property: 



Hjpq 



ipqtj 



Hj qp, 



(2.20) 



which allows us to identify M with a symmetric linear transformation on the space of symmetric 
d x d matrices. It also satisfies the positivity property (positive or negative definiteness) on the 
space of symmetric matrices [H [11] . 

Let U be the background solution associated with (Ao, /io, Po) in ^, i.e., 



(£;Wo + Pow 2 )U = 0, on O, 
dU 



dv 



= g 



on dfl, 



(2.21) 



G 



Then, the following result can be obtained using analogous arguments as in [71 [5] ; see [3] . Here 
and throughout this paper 

d 

A : B = "<A 

for matrices A = (a,ij)fj—i and B = {bij)fj =1 . 

Theorem 2.2. Let u be the solution to (|2.4j) , U be the background solution defined by (|2 . 2 1 1) 
and PqUJ be different from the Neumann eigenvalues of the operator —C\ ^ in {L 2 (Vt)) . Let 
D be given by (|2.1[) and the conditions (I2.2[) and (|2.3|) are satisfied. Then, for u>8 <C 1, i/ie 
following asymptotic expansion holds uniformly for all x € c?Sl : 

u(x)-U(x)=-5 i (vU(z fl ):M(B)V ! „N u (x ) z a ) (2.22) 

+cj 2 (p - pi)|S|N"(x,z a )U(z a )) + 0(<5 d+1 ). 
As a direct consequence of expansion (|2.22[) and Lemma [2TT1 the following result holds. 



Corollary 2.3. Under the assumptions of Theorem 

[u-U](x)= ( S [i (vU(z fl ):M(B)V Zo r^(x-z a ) (2.23) 

W{ Po - pi)|B|rSf(x - z a )U(z a )) + 0(5 d+l ) 
uniformly with respect to x € <9f2. 

Remark 2.4. M^e ftaue made use of the following conventions in (|2.22[) and (|2.23p : 



(VU(Z„) : M(B)V Zo N w (x,z a )) fc = ]T (ftU 3 (z fl ) £ m ijpg d p N£ g (x, z a ) 



and 



(N-(x, z a )U(z a )) fc = ^ N£.(x, z a )U 2 (z a ). 



i=l 



3 Imaging small inclusions using TD 

In this section, we consider a filtered quadratic misfit and introduce a TD based imaging 
functional resulting therefrom and analyze its performance when identifying true location z a of 
the inclusion D. 

For a search point z s , let u z s be the solution to (|2.4[) in the presence of a trial inclusion 
D' = 5'B' + z s with parameters (A^, p[, p\), where B' is chosen a priori and 8' is small. Assume 
that 

dA' 1 +2/x' 1 >0, /ii>0, (Ao-Ai)(^,-Mi)>0. (3.1) 

Consider the elastic counterpart of the filtered quadratic misfit proposed by Ammari et al. in 
[5], that is, the following misfit: 



£/[U](z s ) = i 



Oil 



1 



1 - I u z s _ u mcas ](x) 



dcr(x). (3.2) 
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As shown for Helmholtz equations in [3] , the identification of the exact location of true inclusion 
using the classical quadratic misfit 

£ [U](z s ) = \f |(u zS - u meas )(x)| 2 d<7(x) (3.3) 
1 Jan 

cannot be guaranteed and the post-processing of the data is necessary. We show in the later 
part of this section that exact identification can be achieved using filtered quadratic misfit £f. 

We emphasize that the post-processing compensates for the effects of an imposed Neumann 
boundary condition on the displacement field. 



3.1 Topological derivative of the filtered quadratic misfit 

Analogously to Theorem 12.21 the displacement field u z s , in the presence of the trial inclusion 
at the search location, can be expanded as 



u z s(x) - U(x) = -{5') d VU(z 6 ) : M'(S')V z sN"(x,z 6 ) 



! (p -pi)|S'|N"(x )Z s )U( Z s ) )+0((S') 



(3.4) 



for a small 5' > 0, where M'(B') is the EMT associated with the domain B' and the parameters 
(Ao, i^o'i A'i, /ii). Following the arguments in [5], we obtain, by using Corollary 12 .31 and the jump 
conditions (|2.14[) . that 



S f [U](z s ) = \ [ 
1 Jon 



I-K&) [U-u mcas ](x) 



dcr(x) 



+ (S') d ^e |VU(z s ) : M'(S')Vw(z s ) + uj 2 {p - pi)|S'|U(z s ) • w(z s )} 
+0 ((55') d ) + O ({5') 2d ) , (3.5) 
where the function w is defined in terms of the measured data (U — u mea s) by 



w(x) = St 



-I - K% [u 



U] 



(x), x e ft. 



(3.6) 



The function w corresponds to backpropagating inside f2 the boundary measurements of U 
Umeas- Substituting (|2.23|) in (|3.6j) . we find that 



w(z b ) = S d VU(z ) : M(B) / T%(z b - x)V Za r^(x - z a )d<r(x) 



+u 2 (po-pi)\B\ / r^(x-z a )r^(x-z s )da(x) U(z a ) )+0(6 d+1 ). (3.7) 
l Jdn J ' 

Definition 3.1. (Topological derivative of£f) The TD imaging functional associated with £f 
at a search point z £ fi is defined by 



Xtd[U](z s ) := - 



d£f[V}( z S ) 
d(5>) 



Ad 



(3.8) 
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The functional Ztd [U] (z s ) at every search point z s £ fl synthesizes the sensitivity of 
the misfit £j relative to the insertion of an elastic inclusion D' = z s + S'B' at that point. 
The maximum of Ztd [U] (z s ) corresponds to the point at which the insertion of an inclusion 
centered at that point maximally decreases the misfit £f. The location of the maximum of 
Ztd [U] (z s ) is, therefore, a good estimate of the location z a of the true inclusion, D, that 
determines the measured field u moas . Note that from (|3.5p it follows that 

Z TD [U](z s ) = -3?e{vU(z s ) : M'(B')Vw(z s ) + uj 2 {p a - p' x )\B'\\J(z 3 ) • w(z s )}, (3.9) 

where w is given by (|3.7[) . 



3.2 Sensitivity analysis of TD 

In this section, wc explain why TD imaging functional Ztd may not attain its maximum at the 
location z a of the true inclusion. Notice that the functional Ztd consists of two terms: a density 
contrast term and an elasticity contrast term with background material. For simplicity and for 
purely analysis sake, we consider two special cases when we have only the density contrast or 
the elasticity contrast with reference medium. 

3.2.1 Case I: Density contrast 

Suppose Aq = Ai and po = fix- In this case, the wave function w satisfies 



w(z s )~£ d W 2 (p -Pi)|i?| 



on 



r^(x-z Q )r^(x-z s )da(x) 



U(z a ) 



Consequently, the imaging functional Ztd at z s e fl reduces to 



Z TD [U] (z 6 ') =s Cw 4 3f?e(u(z i '') 



on 



r£(x-z a )r£(x-z s )Ar(x) U(z a 



where 



(3.10) 

(3.11) 
(3.12) 



C = 5 d (p -p' 1 )(p - Pl )\B'\\B\. 
Throughout this paper we assume that 

(Po - p[){po - Pi) > 0. 

Let us recall the following estimates from [31 Proposition 2.5], which hold as the distance 
between the points z s and z a and the boundary <9f2 goes to infinity. 

Lemma 3.2. (Helmholtz - Kirchhoff identities) For z s ,z a £ far from the boundary dil, 
compared to the wavelength of the wave impinging upon £1, we have 



[ r- Q (x-z a )r^ Q (x-z s )d ( T(x) ~ - — 3m{r- Q (z s -z a )} 
Jon c a w 

f ^(x-z^r^fx-z 5 )^) ~ o, a ±p. 
J on 



9 



Therefore, by virtue of (|2.10[) and Lemma 13.21 we can easily get 

1 , o ,1 



I TD [U] (z a ) - -Ci/Se U(z 6 ) 



3m <j — r£ P (z s - z a ) + — r£ s (z s - z a ) }> U(z a ) 

c P cs 



3.13) 



Let (ee l; eg 2 , . . . , e$ n ) be n uniformly distributed directions over the unit disk or sphere, 

(3.14) 



and denote by Jjf and Uf respectively the plane P— and S— waves, that is, 



3 3 

Uf(x) = e lKpx - ee i e 9] and U?(x) = e 



for d = 2. In three dimensions, we set 

U| ; (x) = e ^ x - e ",e^, J = 1,2, 

where (eg . , e e 1 , e^.' 2 ) is an orthonormal basis of R 3 . For ease of notation, in three dimensions, 

I TD [Uf](z s ) denotes Eli ^td[U| ; ](z s ). 
We have 

1 7T . „ 

(3.15) 



-V <" x ' ' ~ -4(— ) d - 2 3mG^(x) 



for large n; see, for instance, [5J. The following proposition holds. 



Proposition 3.3. Let U™ &e defined in (|3.14p . where j = 1, 2, ■ • • , n, /or n sufficiently large. 
Then, for all z s £l! /ar /rom 9f2, 



1 _IL 

I^I TD [Uf](z s )^4 M0 C W 3 (^) d - 2 (^) 2 
n t— ' up Kp 



1 

cp 



3m{r^ P (z s -z a )}| 



{r- p (z s -z a )} :3m{r^ s (z s -z Q )} 



,(3.16) 



and 



1 ™ 1 

- ]>> D [Uj](z S ) c 4 M0 C W 3 (-) d - 2 - |SJm {r- s (z s - z Q )} 



+— Sm {r^ P (z s -Za)}: 3m {r£ 5 (z s - z Q )} 



,(3.17) 



where C is given by (|3.12[) . 
Proof. From (|3.15|) it follows that 



1 ™ 



e ,KpMl >e fi) ®e fij ~ 4( — j'^Sm ( ^D x G P (x) 

3=1 v ^ 



-4 W (f)"(^) 2 Sm{r^(x)} 



(3.18) 
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and 



1 n 1 n 

\ A ksx-efl. _J_ „ J_ * \ A iKqx-efi. / T „ \ 



3=1 



3=1 



4 (_ )^9m I 2 + _D X )G%(x) 



-4 w (-) d - 2 3m{r^(x)} 



(3.19) 



where the last equality comes from (|2.11[) . Note that, in three dimensions, (|3.19p is to be 
understood as follows: 



n 2 



(3.20) 



3 = 1 i=l 



Then, using the definition of U^ 3 we compute imaging functional Itd for n plane P— waves as 

i " i ™ r r 

-Vl TD [Uf](z s ) = C W 4 -y>eUf(z s ). / r^(x-z a )r^(x-z s )d t T(x)Uf(z a ; 



n 



in p (z s — z a ) -en . 



3=1 



{-r^, P (z s -z a ) 

I. Cp 

-T^ s (z s z a )W 



1 ™ 



3=1 



3m <j — r£ P (z s - z a ) + — r£ s (z s - z a ) 

Cp ' cs 



Here we used the fact that eg j ■ Aeg j = eg j <£> eg j : A for a matrix A, which is easy to check. 
Finally, exploiting the approximation (|3 . 18[) . we conclude that 



1 



« Z — ' Kp ftp 



3=1 



-^|3m{r- p (z s -z a )}| 2 



-3m {r£ P (z s - z a )} : 3m {l^ s (z s - z a )} 



Similarly, we can compute the imaging functional Ztd for n plane 5— waves exploiting the 
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approximation (|3.19p , as 

i n 1 n r r 

-Vl TD [Uf](z s ) = 0, 4 -V5ieUf(z s ). / TY(x- Za )r^(x-z s )da(x)Uf(z a ) 



3=1 



3=1 



n 

-Cw 3 -V 5Re 



zk s (z —z a )-e@. _L 



3=1 



|-r- P (z s -z a ) 

t Cp 

-r- s (z s -z Q )}e^ 



-Cw 3 3ie 



3 = 1 



1 ™ 

.{— r- P (z s - Za ) + -r^ s (z s - Za )| 



4 Mo Cu, 3 (— ) c 



|3m{r^ s (z s -z a )}| 



__ 5m {r^p(z s - z Q )} : 3m {r£ 5 (z s - z Q )} 



In dimension 3, one should use (|3.20[) to get the desired result. This completes the proof. □ 
From Proposition 13 . 31 it is not clear that the imaging functional 2td attains its maximum 

^ n 1 71 

at z a . Moreover, for both — \^2rD[Uf ](z s ) and — V^2td[UT ](z s ) the resolution at z a is 

3=1 3 = 1 

not fine enough due to the presence of the term 3m {Tq P (z s — z a )} : 3m s (z s — z a )} . 

_^ n 1 n 

One way to cancel out this term is to combine — y^lTD[U"f](z s ') and — Y^2rD[Uf ](z s ) as 



follows: 



3=1 



- E ( CS (-) d - 2 (-) 2 ^TD[Uf ](«*) - C p(-^)^ 2 I TD [Uf](z 

3=1 

However, one arrives at 



3=1 



1 



3=1 



- E c s (-^) d - 2 (-^) 2 I TD [Uf }(z s ) - cp(^-) rf - 2 I TD [Uf](z 5 ) 

7T Kg 7T J 



o u » f |3m {r- P (z s - Za) }\ 2 -2L | 3m {r- s (z s - z a )}| 2 Y 



which is not a sum of positive terms and then can not guarantee that the maximum of the 
obtained imaging functional is at the location of the inclusion. 
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3.2.2 Case II: Elasticity contrast 

Suppose po = pi. Further, we assume for simplicity that M = M'(B') = ~M(B). From Lemma 
13.21 we have 

f V Za Tf(x-z a )V zS r^(x-z 5 )da(x) ~ -3m{V Za V z3 r^(z s -z a )} 

Jan csu (32i) 

— L Sto {v Za v zS r^ P (z s - Za )} . 

Then, using (j3~7|) and (|3T2"T]) . I T d [U] (z s ) at z s € ft becomes 
Itd[U](z 5 ) = -5 d $teW(z s ) : MVw(z s ) 

^ / V Za ry(x-z Q )V zS r^(x-z 5 )dcr(x) :MVU(z a ) 
Jon 

M V 2 (^m{f${z s -<)}) : MVU(z Q ) 



= J d 3?eVU(z- 
~ — !ReVU(z 



where 



r-(z 5 - z a ) = — r^ P (z 5 - Za ) + — r^ s (z 5 - z a ). (3.23) 

Cp C5 ' 



Let us define 

J Q ^(z s ) := (lS m [(V 2 r^ Q )(z s -z a )]) : (lSm[(V 2 ry(z s -z„)]) T , (3.24) 

where A T = (Af-nj) if A is the 4-tensor given by A = (A^;). Here A : B = J^ijki ^ijkiBijki for 
any 4-tensors A = {A ijkl ) and B = {B ijM ). 
The following result holds. 

Proposition 3.4. Lei U" be defined in (|3.14p . where j = 1, 2, ■ • • , n, /or n sufficiently large. 
Let J a> p be defined by (|3.24p . Then, for all z s G SI far from dft, 



V> D [Uf ](z s ) ~ 4<5^(^) d - 2 (^) 2 (^J P ,p(z s ) + - J s ,p(z s )) (3.25) 

* — ' W Kp Kp \Cp '"o 



and 



n ^— J cj k s v c s c P 

Proof. Let us compute Itd for n plane P— waves, i.e. 



1 ™ 11 

* y> D [Uf](z s ) ^ 4«5^(^) d - 2 (^J s , s (z s ) + -J S p(z s )). (3.26) 



1 n 1 

-^I TD [Uf](z s ) = --Ke^VUf(z s ) :M[3m{(V 2 ^)(z s -z a )} :MVUf(z a ) 



n ' — ' J u) n 

3 = 1 3=1 



i _ 



Cp 71 



3=1 



,{v 2 r^(z s -z a )} :Me 9j <g»e e ,) . (3.27) 
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Equivalently, 

1 n 1 n d d 

i$> D [Uf](z S ) = S^-Sfte^e^-^ £ £ 

i,k,l,m— 1 i' ,k' ,1' ,m' = 1 

$,fff) (z s - z a )) } m iWi , fc , ( 3 - 28 ) 



n * — ' J c P n 

3=1 3 — 1 i,k,l,m—l i ,k ,1 ,m — 1 



mk 

where the matrix A 6j = (A^)^ is defined as A 9j := e#. ® e#.. It follows that 

1 n d d 

-^I TD [Uf](z s ) = S d $te Yl £ m imlfc m /w ^^m[((92,f^)(z s -z a ) 

j — 1 i.kd,m—l i r ,k' ,1' .m' — l 

' (V 1 



mk' 



X^Hf^ Za) ^4i4L>)- (3-29) 

Recall that for n sufficiently large, we have from (13.181) 

i£ e *K PX .e 8 „ ^(JLf-^HlfZmiTZ , P (x)} 

n * — ' Kp ftp 

(with the version p.20[) in dimension 3). Taking the Hessian of the previous approximation 
leads to 

1 n 2 

lv- e « P x. ee „ 4jUo £P(JL)d-2 ( ^£ ) 2 3m { V 2 r ^ (x) } 

n * — ' uj z Kp Kp 

3=1 

c 4 M0 ^V( — )"3m.{V 2 r^(x)}. (3.30) 

W Cg Kp 

Then, by virtue of (|3.18[) and (J3T30J) , we obtain 
iyi TD [Uf](z s ) ~ s*^—)*- 2 ^ 2 

n UJ Kp Kp 

j — l i,k,l,m—li',k',l',m' = l 



,! 



Vx T D[Uf](z s ) ~ 6 d ^(—) d - 2 (^f V V mi mik m Vm>vk , 

— ^ t,) K r> ftp ' ^ ' 

i,k,l,m— 1 i' ,k' ,1' ,m 

3™ { ((^fft) (z s - z a )) } 3m {((df^p) (z s z a )) m , k } 



mk' 

s -^/- 2 2 t (e ■*-*»{((**»)(-•--.))«}) 

i,k,i f ,k f — 1 \ /,m— 1 / 



X 



\ i',m'=l 

Therefore, by the definition (|3.24j) of J Qj /3, we conclude that 



£ m^ ia ^m{((^r^ P )( Z s - z a )) m , k } 



iVl TD [Uf](z s ) ~ ^(^) d - 2 (^) 2 (M9mfv 2 f^(z s -z a )|) 

71 ' J UJ Kp Kp \ l J / 



M3ro{V 2 r^ P (z s - z a )} 



7' 



^^(-) d ~ 2 (^) 2 ( - JpA* s ) + -JsA* s ) 

UJ Kp Kp \Cp C$ 
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Similarly, consider the case of plane S— waves and compute Xtd for n directions. We have 
-^2T TD [Uf](z s ) = --& e ^VUf(z s ) : M (Qm { (V 2 r£) (z s - z a )} : MVUf (z a ) 



n e — ' w n 

3=1 3=1 



£> 3 = 1 



3m{(V 2 r£)(z s -z a )} iMej-Je^ 



1 n d d 

- ^>E^ s(zS ~ z ^ E E 

^ j—l i,k } l,m=l i',k' ,1' .in' — 1 

xSmjf^f^tz^zjj^Jm^wB^ (3.31) 
where the matrix T$ 0j = (B^ )^ is defined as B 0j = ® ej- . It follows that 



1 n d d 

-E^TD[Uf](z s ) = ^ E E ">, mi)t m w «-3mfe (f^(z 5 -z Q )) 



j — 1 i,k,l,m— 1 i' ,k f ,1' ,m' — 1 



^l^^-^BtiB^j. (3.32) 
Now, recall from (|3.19[) that for n sufficiently large, we have 

iE^ sx - e "H ® e i - -Vo(-) d - 2 Sm{r- s (x)} . 

3=1 

Taking the Hessian of this approximation leads to 

i n 2 

£ W« S x-e 9 ^ I ®e ej ®et =-4 Mo %— ) d - 2 3m{V 2 r^(x)} , (3.33) 

3 = 1 

where we have made use of the convention 

(V a rfo) tfH <UH^) jr 
Then, by using (|3.19[) . (|3.33[) and the similar arguments as in the case of P— waves, we arrive 
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at 



( ) > J > J mimik mi' m'i'k' 

U us , ,, 

i,k : l,m—l i ,k ,1 ,m — 1 

x3m{((9f^)) mfc/ (z s -z a )} 

x3m{((^r^)) m , fc (z s -z a )} 

" ) d " 2 (M9m{(V 2 f^)(z s - Za )}): 

(M^{(v 2 r- s )(z s - Za )}) T 

* 5^(^)^(J-Jp, s ^) + J-J SiS ^)). 

This completes the proof. □ 

As observed in Section 13.2.11 Proposition 13.41 shows that the resolution of Itd deteriorates 
due to the presence of the coupling term 

J PjS (z 5 ) = (MSm{(V 2 r^ s )(z s -z a )}) : (lS m {(V 2 r^)(z s -z a )}) T (3.34) 
3.2.3 Summary 

To conclude, we summarize the results of this section below. 

- Propositions 13.31 and 13.41 indicate that the imaging function Itd may not attain its max- 
imum at the true location, z Q , of the inclusion D. 

- In both cases, the resolution of the localization of elastic anomaly D degenerates due 
to the presence of the coupling terms 3m {Tq P (z s — z a )} : 3to {Tq s (z s — z a )} and 
Jp i s(z s '), respectively 

- In order to enhance imaging resolution to its optimum and insure that the imaging func- 
tional attains its maximum only at the location of the inclusion, one must eradicate the 
coupling terms. 



^]>> D [Uf](z s ) * 



4 Modified imaging framework 

In this section, in order to achieve a better localization and resolution properties, we introduce a 
modified imaging framework based on a weighted Hclmholtz decomposition of the TD imaging 
functional. We will show that the modified framework leads to both a better localization (in the 
sense that the modified imaging functional attains its maximum at the location of the inclusion) 
and a better resolution than the classical TD based sensitivity framework. It is worthwhile 
mentioning that the classical framework performs quite well for the case of Hclmholtz equation 
[5] and the resolution and localization deteriorations are purely dependent on the elastic nature 
of the problem, that is, due to the coupling of pressure and shear waves propagating with 
different wave speeds and polarization directions. 
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It should be noted that in the case of a density contrast only, the modified imaging func- 
tional is still a topological derivative based one, i.e., obtained as the topological derivative of 
a discrepancy functional. This holds because of the nonconversion of waves (from shear to 
comprcssional and vice versa) in the presence of only a small inclusion with a contrast den- 
sity. However, in the presence of a small inclusion with different Lame coefficients with the 
background medium, there is a mode conversion; see, for instance, [21]. As a consequence, the 
modified functional proposed here can not be written in such a case as the topological derivative 
of a discrepancy functional. It is rather a Kirchhoff-type imaging functional. 

4.1 Weighted imaging functional 

Following [3J, we introduce a weighted topological derivative imaging functional Z\y, and justify 
that it provides a better localization of the inclusion D than Ztd ■ This new functional I\y can 
be seen as a correction based on a weighted Helmholtz decomposition of Ztd- I n fact, using the 
standard L 2 -theory of the Helmholtz decomposition (see, for instance, [16]), we find that in the 
search domain the pressure and the shear components of w, defined by (|3.6p . can be written as 

w = V x Vw + V0 W . (4.1) 

We define respectively the Helmholtz decomposition operators H, F and Ti by 

U P [w] := V0 W and U s [w] := V x Vw (4.2) 

Actually, the decomposition w = 7i p [w] + H s [w] can be found by solving a Neumann prob- 
lem in the search domain [16] . Then we multiply the components of w with cp and cs, the 
background pressure and the shear wave speeds respectively. Finally, we define I\y by 

Zw[U] = cpSRe j-VH p [U] : M' {B')VH p [w} + lj 2 - lj \B'\H P [U] -H p [w]| 

+ c s Ue |-VH S [U] :M'(B')\/n s [w}+uj 2 Oj^-l) |B'|H S [U] ■ H 5 [w] | .(4.3) 

We rigorously explain in the next section why this new functional should be better than imaging 
functional Itd • 

4.2 Sensitivity analysis of weighted imaging functional 

In this section, we explain why imaging functional 2Av attains its maximum at the location z a 
of the true inclusion with a better resolution than Itd- In fact, as shown in the later part of 
this section, I\y behaves like the square of the imaginary part of a pressure or a shear Green 
function depending upon the incident wave. Consequently, it provides a resolution of the order 
of half a wavelength. For simplicity, we once again consider special cases of only density contrast 
and only elasticity contrast. 

4.2.1 Case I: Density contrast 

Suppose Ao = Ai and fiQ — /ii. Recall that in this case, the wave function w is given by (13.101) . 
Note that H 01 ^] = T% al a G {P, S}. Therefore, the imaging functional lw at z s e £1 turns 
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out to be 

l w [XJ](z s ) = Cw 4 5ie f c P H p [U](z s ) 



T-(x - z a )r^ P (x - z s )da(x) U(z Q ) 



c s n s [v](z s )- / r-(x-z )re s (x-z s )<fa(x) u(z a ) . (4.4) 



By using Lemma 13721 we can easily get 

I w [U](z s ) ~ -Cu 3 Me\H p [V](z s )- [3m{q p (z s -z a )}U(z„) 

+H s [U](z s ) • fern {r^ s (z s - z a )} U(z a ) 



(4.5) 



Consider n uniformly distributed directions (e$ 1 , eg 2 , . . . , eg n ) on the unit disk or sphere for 
n sufficiently large. Then, the following proposition holds. 

Proposition 4.1. Let U" be defined in (|3.14|) . where j = 1, 2, ■ • • , n, for n sufficiently large. 
Then, for all z s E fl far from dft, 



1 



- $>w[Uf }(z s ) c ^C^i^r-^f |3m {T" 0P (z s - z a )}| 2 , 

« J Kp Kp ' ' 



(4.6) 



3=1 



and 



1 n 

- ^I w [Uf](z s ) ~ 4 Mo O, 3 ( — ) d ~ 2 |3m {T^ s (z s - z a )}| 2 , (4.7) 



n <r-. ' k s 

3 = 1 



where C is given by (|3.12[) . 

Proof. By using similar arguments as in Proposition 13.31 and (|4.5[) . we show that the weighted 
imaging functional Xw for plane P— waves is given by 



1 



]Tl W [Uf](z s ) = -C 3 --3?e^Uf(z s ). 3m{r^(z s -z„)}Uf(z n ) 



3=1 



3=1 



1 " 



For n plane S 1 — waves 

n 

-5>w[Uf](z S ) 



3 = 1 



3=1 

4 / , C^ 3 ( — )^ 2 (--i) 2 |3m {r" P (z s - z a )}| 



n 

' CW3 - E ' [ 5 ™ - Za)} Uf (Z a )] 

3 = 1 



1 



3 = 1 

n 



{T^ s (z s - z a )} ei. 



4 Mo C W 3 ( — ) d ~ 2 \<&m {r^ s (z s - z a )}f , 

K S 
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where one should use the version (|3.20[) in dimension 3. 



□ 



Proposition 14. 1 1 shows that Zw, attains its maximum at z a (see Figure [T]) and the coupling 
term 3m {Tq p (z s — z a )} : {Fq s (z s — z a )}, responsible for the decreased resolution in 
Ztd, is absent. Moreover, the resolution using weighted imaging functional Zw is the Raylcigh 
one, that is, restricted by the diffraction limit of half a wavelength of the wave impinging upon 
O, thanks to the term |3t71 {Tp Q (z s — z a )}| . Finally, it is worth mentioning that Zw is a 
topological derivative based imaging functional. In fact, it is the topological derivative of the 
discrepancy functional cs£f [U s ] +cp£f [U p ] , where XJ S is an S'-plane wave and U p is a P-plane 
wave. 





Figure 1: Typical plots of |3m{r£ s (z s - z a )}\ 2 (on the left) and \%m{T% P (z s - z a )}\ 2 (on 
the right) for z a = and cp/cs = VTl. 



4.2.2 Case II: Elasticity contrast 

Suppose po = pi and assume for simplicity that 
imaging functional Zw reduces to 



'(£?') = M(B). Then, the weighted 



Z w (z s ) 



c P VH r [U{z b )] : MVH^[w(z i )] + c s VH b [U {z b )} : MV^w^)] 



c P V"H p [U(z s )] : M( / V Za r^(x-z a )V zS r^ P (x-z s )dcr(x) : MVU(z Q ; 

an 



S VH s [U(z s )] : Ml / V Zn r^(x - z a )V zS r^ jS (x - z s )d<r(x) : MVU(z a ) 



-6 d 



Wn p [U(z s )} :M(Sm{(V 2 ry( Z s -z a )} :MVU(z„) 



+VM s [U(z s )]:M3m{(V 2 r^)(z s - Zo )} : MVU(z a ) 



(4.8) 



We observed in Section 13.2.21 that the resolution of Ztd is compromised because of the 
coupling term Js,p(z s ). We can cancel out this term by using the weighted imaging functional 
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Iw- For example, using analogous arguments as in Proposition 13.41 we can easily prove the 
following result. 

Proposition 4.2. Let U" be defined in (|3.14p . where j = 1, 2, • • • , n, for n sufficiently large. 
Let J a _p be defined by Q3.24p . Then, for all z s £ far from dtt, 



i^Iw[U«](z s )^ 4( 5^(^) d - 2 (^)V Q . Q ( Z s ), ae{P,S}. 



3=1 



(4.9) 



It can be established that Iw attains its maximum at z s = z a . Consider, for example, the 
canonical case of a circular or spherical inclusion. The following propositions hold. 

Proposition 4.3. Let D be a disk or a sphere. Then for all search points z s G f2, 
Jp P (z s ) = a 2 V 2 (3mr^ P )(z s -z a ) 2 + 2a6 A(3mr^ P )(z s -z Q ) 



+b 2 ATr(3mry(z s -z a ) , 

where Tr represents the trace operator and the constants a and b are defined in (12.191) . 
Proof. Since 



fv 2 r" 



it follows from (|2~T9|) that 



(MV 2 r^ 



' ijkl 



iij pq (V 2 TQ ,p) pqkl 



(4.10) 



(4.11) 



(4.12) 



3=1 

= \^ ((VT^e,^ + (VT8>e,)y) + 69 fc V • ((r^e,))^, (4.13) 

where e; is the unit vector in the direction X;. 

Now, since Tq p e; is a P— wave, its rotational part vanishes and the gradient is symmetric, 
i.e., 

Vx(M = and (Vr^e,) . . = (Vr^) .. = (Vr^ F e^. • (4.14) 
Consequently, 

W- ((r£pe,)) = V x (v x (r^pe ; )) +A(r-pe ; ) =A(r- pei ), (4.15) 
which, together with (|4. 13[) and (|4.14[) . implies 

MV 2 r^.p =aV 2 r^ P + 6I 2 ® Ar^p. (4.16) 
Moreover, by the definition of T% p , its Hessian, V 2 Tq p , is also symmetric. Indeed, 

(V 2 r- P ) T = d ki {Tl P ) = -^d ktjl G P = h 2 Kp) (4.17) 



ijkl 
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Therefore, by virtue of (|4.16[) and (|4.17[) , Jp,p can be rewritten as 

JpA* S ) = (aQm{(V 2 r^ P )(z s -z a )}+6I 2 ®3m{(Ar^ P )(z s -z a )} 
: (a3m{(V 2 r^ P )(z s - z a )} + fe3m{(Ar^ P )(z s - z a )} ® I 
Finally, we observe that 

V 2 3m{r^ P }) : (v 2 9m{r^ P }) T = V 2 9m {i^p} 



(4.18) 



(4.19) 



V 2 3?7i{r^p} : (l 2 ® A3m{r£ P }) = V 2 Zm{T% P } : (ASm{r^ P } ® I 2 ) 

d 

= £ (Sm(4r^)J%A3m(r^) fci 

fc,i=l \ 4=1 / 



E ( A ^( r 0,p) fei ; 
k,l=l 

A3m{ry 2 , 



(4.20) 



and 



(l 2 ® A9fm{rSf (P }) : (AQfm^p} ® I 2 ) = E % A9 m (r^) H 4,A3m (rg^ 

i,j,k,l—l 



a 

= E A3m ( r 0,p)M AS ™( r 0,p), 

i,fc=l 



ATr(9m{r^ P }) i 

We arrive at the conclusion by substituting (|4. 19[) . (|4.20p and (|4.21j) in (|4.18|) . 
Proposition 4.4. Let D be a disk or a sphere. Then, for all search points z s G £7, 



(4.21) 
□ 



Mo 



h 's 



V 4 3m {G£(z s - z a )} J" + |v 2 3m {G£(z s - z Q )} 



-^b m {^(z' 



»)} 



Mo 



Kg 



4r E l^ 5 ™ 

kl,k^l 

{G£(z s -z a )} 



{G%(z s - Za )}| 2 + i^|v 2 Sm{G^ 



(4.22) 
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where a is the constant as in (|2.19[) . 
Proof. As before, we have 

( MV 2 r^ s 



ijkl 



dik (rfo)., ■ <>ji- (T% s ) a )+bd k V- ( (rgf iS e,) )S. 



and 



(4.23) 



(4.24) 



Here we have used the facts that s ei is a 5— wave and, Tq s and its Hessian are symmetric, 

d ik (T^s) j{ = d M (i^) = d ik (Tls) x . = d ki (r- s ) ( . . (4.25) 
Substituting, P~2"5")) and P~2"4]l in (f3T2"4"|) . we obtain 

2 d 

^s(z s ) - T E 3m 



i,j,k,l— 1 



(^r^ s )(z s - + ((3 Jfe r^)(z s - Za )) i; } 
x3m{((&ry(z s - Za )) j7 + ((^r^ s )(z s - z a )) .J 

^(T 1 (z s ) + 2T 2 (z s ) + r 3 (z S )), (4.26) 



where 



T!(Z S ) 



r 2 (z s ) 



(z s - 



r 3 (z S ) 
Notice that 



£ (sm{(&r^) 3! (z s - z Q )}) (9fm{(a ifc rg, s ) .,( 

z ,j,k,l = 1 

X; (s™{(a 4fc r^ s ) j; (z s -z Q )}) (sm{(^r- s ). fe ( 

d 

]T (s>m{(a,- fc rsf iS ) a (a s -a tt )}) (3m{(a,r- s ) jfc (z s -z a )}) 



i,j,k,l— 1 



•)}) 
«)}) 



i 



3m{r^ s (x)} = 5.(412 +© x )3fm{GS(x)} ) 



and Sm{Gg} satisfies 

A3m {G5} (z 3 - z a ) + k|3to {G£} (z 5 - z a ) = for z s ^ z a . 
Therefore, the first term T\ can be computed as follows 

Ti(z S ) = |v 2 (3mr^ s )(z s -z a 

d 



(4.27) 



g (»m Gg) (z s - z a )) 2 + 4Sjt (d ik (9fm Gg) (z 5 - z a )) 

+24^9^ (3m Gg) (z s - z„)fi|y M (3mG^) ( 



i^j,k,l—l 



(z s - Z„ 
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We also have 



53 ^ildik (S*m Gg) (z s - z a ) (d ljkl (3m Gg) (z s - z a ; 

i,j,/c,Z— 1 



d ( d \ 

= 2 ^ (ft fe (3m Gg) (z s -z a )) 5 ife ^ S H (3m Gg) (z 5 - z a )\ 

i,k=l V Z=l / 

d 2 

= -2«| ^ (9, t (3mG^)(z s -z a ))", 



i,k=l 



and 



53 S iX (d lk (3m Gg) (z 5 - Zo )) = d 53 (a ifc (3m Gg) (z s - z a ; 
Consequently, we have 



i,fc=l 



T!(z s ) = V 2 (3mr- 5 )(z s -z Q ) =^V 4 (3mG^)(z s -z„ 



■^E- fc=1 (5 lfe (^Gg)(z s -z a ) 



Mo 

Estimation of the term T2 is quite similar. Indeed, 

d 



(4.28) 



i.j : k : l — l _ 



9, jH (3mG^)(z s -z a ) 

+24^,^ (3m Gg) (z 5 - z a )% H (3m Gg) (z 5 - z a ) 
+4^ fc (d 4fc (3mGg)(z 5 - z a )) (a i; (3mGg)( 



z s - z„ 



Finally, using 

d 



53 ^^ fc (a ifc (3mGg)(z s -z Q ))(a 4i (3mGg)(z s -z a )) = 53 (d 4fc (3mGg)( 



z s - z„ 



i,j,k,l— 1 

we obtain that 



;,fe=i 



T 2 (z*) = — 



2 1 



V 4 (3mGg)(z s - z a ) - — V 2 (3mGg)(z s - z a ) 



(4.29) 



Similarly, 



9 yfc; (3mGg)(z s -z a ) 



+24^9i fe (3m Gg) (z s - z a ) (d ijkl (3m Gg) (z 5 - z a ) 



4*«*;fc (3m Gg) (z 5 -*„))( d a (3m Gg) (z s - z a ) 
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By virtue of 



tufa fe(3mGs)(z S - *„)) (d a {^mG%)(z s - z a )) 

i,j,k,l=l 



a 

= ( 9kk G S) - Z -)) { 9 H G S) - Za)) 

i,fc=l 



wc have 



Mo 



2 2 



V 4 (3mG^)(z s - z a ) - — V 2 (3mG^)(z s - z a ) 



3mG^(z s -z a ) 



(4.30) 



We conclude the proof by substituting (|4.28|) . (|4.29[) and (|4.30[) in (|4.26l) and using again 
(OTP . □ 

Figure [5] shows typical plots of J ayCc for a S {P, 5 1 }. 




-4 -2 2 4 6 



Figure 2: Typical plots of J55 (on the right) and Jpp (on the right) for z a = and cp/cs 

VTT. 



5 Statistical stability with measurement noise 

Let \5j and be as before. Let {Uj} be plane waves. Define 

n 

Iwf[{U,}](z s ) = -^I w [U,](z s ). 

n 3=1 



(5.1) 



In the previous section, we have analyzed the resolution of the imaging functional Twf in the 
ideal situation where the measurement u moas is accurate. Here, we analyze how the result will 
be modified when the measurement is corrupted by noise. 
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5.1 Measurement noise model 



We consider the simplest model for the measurement noise. Let u truG be the accurate value of 
the clastic displacement field. The measurement u mcas is then 

U m cas(x) = U t rue(x) + l'noise(x), (5.2) 

that is the accurate value corrupted by measurement noise modeled as f no i S e(x), x € 90. Note 
that f no ise(x) is valued in C d , d = 2,3. 

Let E denote the expectation with respect to the statistics of the measurement noise. We 
assume that {i>' no i so (x), x G dQ} is mean zero circular Gaussian and satisfies 



E[^ n0 isc(y) ® ^noi SO (y')] = <nsA(y')l2- ( 5 - 3 ) 

This means that firstly the measurement noises at different locations on the boundary are 
uncorrclatcd; secondly, different components of the measurement noise are uncorrelated, and 
thirdly the real and imaginary parts are uncorrelated. Finally, the noise has variance c^ oisc . 

In the imaging functional 2"wf, the elastic medium is probed by multiple plane waves with 
different propagating directions, and consequently multiple measurements are obtained at the 
boundary accordingly. We assume that two measurements corresponding to two different plane 
wave propagations are uncorrelated. Therefore, it holds that 



EK„ isc (y) ® ^oisc(y')] = <o is JjiS y (y')h, (5.4) 

where j and I are labels for the measurements and Sji is the Kronecker symbol. 

5.2 Propagation of measurement noise in the back-propagation step 

The measurement noise affects the topological derivative based imaging functional through the 
back-propagation step which builds the function w in (|3.6p . Due to the noise, we have 



w(x) = S£ 



— ) [U — Utruc — ^noisc] 



(x) = Wtrue(x) + W noisc (x), (5.5) 



for x £ O. Here, w trU e is the result of back-propagating only the accurate data while w no i SO is 
that of back-propagating the measurement noise. In particular, 



w noise( x ) — — Sq 



(x), xe!]. (5.6) 



To analyze the statistics of w no i so , we proceed in two steps. First define 

^noisc,l( x ) = Q-T - ^nj [^noisc](x), x e dQ. (5.7) 

Then, due to linearity, f no isc,i is also a mean-zero circular Gaussian random process. Its 
covariance function can be calculated as 

1 r , _ 1 



E[l/noise,l(y) ® ^noiscl(y')] = ^E[z/ noiso (y) (g) U noiae (y')] - -E[/Cq [f n oise] (y) (g) fnoisc(y')] 



-^E[iy noise (y) <g> fC%[v noiae \(y')] + E[/C£[i/ noise ](y) <g> tC%[v noise ](y>)}. 
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The terms on the right-hand side can be evaluated using the statistics of f no isc and the explicit 
expression of JCq. Let us calculate the last term. It has the expression 



E 



an Jan 



dT; 



dv- 



-(y - x)iv noiso (x) 



dT 



dv- 



^(y'-x'Koise(x') 



da(x)dcr(x ) 



Using the coordinate representations and the summation convention, we can calculate the jfcth 
element of this matrix by 



an Jan 



dv* 
dT 



(y-x) 



Oil 



dv* 
&T% 

a dv x 



L (y-x] 



dT, 
dv x , 

dT 



My'-x') 
°(y'-x 



dv 



EKoiseW ® v noisc (x')}i s d(7(x)da(-x') 



da(x) 



dn 



(y-x)^°(x-y')d<7(x). 
dv* 



In the last step, we used the reciprocity relation 

r£(y-x) = [r£(x-y)] T , 



(5.8) 



for any x, y € K d . 

The other terms in the covariance function of u no ise,i can be similarly calculated. Conse- 
quently, we have 



E[fnoise,i(y) ® ^noisc.l (y')l 



2 2 
_ "noise r /',,'\T "noise 



on 



on 

dv* 



^(y-y') + f iL (y-y') 

dVy> OVy 

ffpu 

x)-^(x-y')da(x). 



(5.9) 



From the expression of Iwf and %, we see that only the Helmholtz decomposition of w meas , 
that is H p [w] and H [w], are used in the imaging functional. Define w Q = W a [w], a <G {P, S}. 
Using the decomposition in (|5.5[) . we can similarly define w" ruo and w^ oiso . In particular, we 
find that 



sc(x) = - / r^ Q (x - y)^noiso,i(y)^( y ), xed. 

JdQ 



This is a mean zero C d -valued circular Gaussian random field with parameters in f2. The jfcth 
element of its covariance function is evaluated by 



EK oi8e (x)®w«. se (x')] jfc = W (ro, M)jl(^, Q (x'-y')) fes E[iwT(y)®^noisca(y') 

l,s 



Is- 



Using the statistics of f no isc,i derived above, we find that 



E[w« oisc (x)®w« oisc (x')] 
Io, a (x-y 



on 



rSU(x-y)rS )Q (y-x')d<7(y) 



dT% , „ . cT 



(y - y') + p-(y - y') r- Q ( y ' - x!)da{y)da{y') 
(any L w y dv y , J 

piru ffru 

r£ Q (x - y)^(y - z)-^(z - y')r^ a (y' - X ')da(z)da(y)da(y>). 
(an) 3 c^z ov z 



2G 



Thanks to the Helmholtz-Kirchhoff identities, the above expression is simplified to 



EK Jx) ® w« oisc (x')] = -^Sm{^ a (x - x')} 



g 2 . f d3m{r" (y-x')i 



ggm{r^(x-y')} _ , 

^ r o,a(y - x )My ) 

9Sm{ry >a (x - z)} 9Sm{r^ a (z - xQ} 
777, 777, da{ ? 
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2c a a; Jqq 

2 n 

^ ^noise / 

(c a uj) 2 J dn dv z dv. 



Assuming that x, x' are far away from the boundary, we have from [3] the asymptotic formula 
that 

dr" (x-y) 

"-"iCcuTZJx-y), (5.10) 



dVy 



where the error is of order o(|x — y| 1//2 d ). Using this asymptotic formula and the Helmholtz- 
Kirchhoff identity (taking the imaginary part of the identity) , we obtain that 



a 2 



E[w° oisc (x) ® w° oi8C (x')] = -^Smfr^fx - x')}. (5.11) 

In conclusion, the random field w I " oise (x),x G Q, is a Gaussian field with mean zero and 
covariance function (|5.1ip . It is a speckle pattern, i.e., a random cloud of hot spots where 
typical diameters are of the order of the wavelength and whose typical amplitudes are of the 
order of cr noi8e /(2^/c a w). 



5.3 Stability analysis 

Now we are ready to analyze the statistical stability of the imaging functional Zwf- As before, 
we consider separate cases where the medium has only density contrast or only clastic contrast. 



5.3.1 Case I: Density contrast 

Using the facts that the plane waves U p 's are irrotational and that the plane waves U s 's are 
solcnoidal, we see that for a searching point z £ O, and a £ {P,S}, 

Z WF [{U«}](z) = c a co 2 - lj |B'|-£>{U?(z) • K tmc (z) + w« noisc (z))}. 

We observe the following: The contribution of {w" truo } are exactly those in Proposition 14.11 
On the other hand, the contribution of {w" noiso } forms a field corrupting the true image. With 
C a := c a uj 2 \B'\(pi/ po — 1), the covariance function of the corrupted image, can be calculated 
as follows. Let z' e 0. We have 

1 " 

Cov(I W f[{U«}](z),Iwf[{U«}](z')) = Cl-^ ]T E[Ke{U? • w^JiefUj" ■ w£ nois J] 

U 3,1=1 

= G l i E ^ {W ■ E Knoi S c(^) ® ^~(^)]UP} • 

3=1 
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To get the second equality, we used the fact that w" noise and w" noise are uncorrected unless 
j = I. Thanks to the statistics (|5.11[) . the covariance of the image is given by 

2 i n 

a i =1 

where e@. = eg . and eg. = e^. . 

Using the same arguments as those in the proof of Proposition l4.il we obtain that 

Cov(X W f[{U^}](z),Z wf [{U5"}](z')) = C' a ^\Zm{r» a (z z')}| 2 , (5.12) 
where the constant 

c" a = c a co^o\Bf(^-in-r-"(^r. 

PO K a K a 

The following remarks hold. Firstly, the perturbation due to noise has small typical values 
of order a no i SC / \/2n and slightly affects the peak of the imaging functional Iwf- Secondly, the 
typical shape of the hot spot in the perturbation due to the noise is exactly of the form of the 
main peak of Iwf obtained in the absence of noise. Thirdly, the use of multiple directional 
plane waves reduces the effect of measurement noise on the image quality. 

From (|5.12[) it follows that the variance of the imaging functional Iwf at the search point 
z is given by 

Var(Z W F[{U«}](z)) = ^^|3m{r£ Q (0)}| 2 . (5.13) 
Define the Signal-to-Noise Ratio (SNR) by 

;[Z W F[{U«}](Za)] 



SNR := 



Vax(l WF [{V?}}(z a ))^> 
where z a is the true location of the inclusion. From (|4.6p . (|4.7[) . and (|5.13[) . we have 



SNR = -¥ |3m{r- a (0)}|. (5.14) 

C'noise 

From (|5.14[) . the SNR is proportional to the contrast |pi — po I an d the volume of the inclusion 
<5 d |.B|, over the standard deviation of the noise, cr no i sc . 

5.3.2 Case II: Elasticity contrast 

In the case of elastic contrast, the imaging functional becomes for z £ Q 

n 

Xwf[{U?}](z) = c a -^VU«(z) : M'(B')(Vw? true (z) + Vw? noise (z)). (5.15) 

11 3 = 1 

Here, w" true and w" noiso are defined in the last section. They correspond to the backpropa- 
gation of pure data and that of the measurement noise. The contribution of w" truo is exactly 
the imaging functional with unperturbed data and it is investigated in Proposition 14.21 The 
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contribution of w" noiso perturbs the true image. For z,z' £ f2, the covariance function of the 
TD noisy image is given by 

Cov(I WF [{U«}] (z) , I WF [{U«}] (z')) 
1 " 

= c '^2 E E[^e{VU«(z) : M'Vw« noisc (z)}3?e{VUr(z') : M'Vw^Jz')}] 

3,1=1 



I n 

=c «^ E ^E[(VU«(z) : M'Vw« noisc (z))(VU«(z') : M'Vw« nois c (z'))] 
1 n 

= c ^E fc { VU "( z ) :M'[E[Vw; oisc (z)Vw; oisc (z')] : M'VUf (z')] } 



in" 



Using (|5.1ip . we find that 



E[Vw^ noisc (z)Vw« noisc (z')] = ^« z V,{r^ a (z - z')}. 

After substituting this term into the expression of the covariance function, we find that it 
becomes 



r a 1 1 

L a u noise 1 



^5Re{vU?(z) : M' 3m{V 2 r^ a (z - z')} : M'VUj(z') } . 



4w In 

i=i 

The sum has exactly the form that was analyzed in the proof of Proposition 13.41 Using similar 
techniques, wc finally obtain that 

Cov(I W f[{U«}](z),Iwf[{U«}](z'))=Mo(-) 3 (— ) d ~\ — ) 2 %^ J a , a (z,z'), (5.16) 

UJ foot f^ot ZiTL 

where J Q Q is defined by (|3.24p . The variance of the TD image can also be obtained from (15.161) . 
As in the case of density contrast, the typical shape of hot spots in the image corrupted by 
noise is the same as the main peak of the true image. Further, the effect of measurement noise 
is reduced by a factor of y/n by using n plane waves. In particular, the SNR of the TD image 
is given by 

5 d y/UQLu , n ns 4\/2n 



SNR= ^ ( — ) 2 ^/j Q , Q (z a ,z a ). (5.17) 

K>a. K>a ^noise 



6 Statistical stability with medium noise 

In the previous section, we demonstrated that the proposed imaging functional using multi- 
directional plane waves is statistically stable with respect to uncorrelated measurement noises. 
Now we investigate the case of medium noise, where the constitutional parameters of the elastic 
medium fluctuate around a constant background. 

6.1 Medium noise model 

For simplicity, we consider a medium that fluctuates in the density parameter only. That is, 

p(x)=p [l+7(x)], (6-1) 
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where p$ is the constant background and po7( x ) is the random fluctuation in the density. Note 
that 7 is real valued. 

Throughout this section, we will call the homogeneous medium with parameters (Ao, po, po) 
the reference medium. The background medium refers to the one without inclusion but with 
density fluctuation. Consequently, the background Neumann problem of elastic waves is no 
longer (|2.2f p . Indeed, that equation corresponds to the reference medium and its solution will 
be denoted by U®. The new background solution is 

(£ Ao , M „+Po^ 2 [l+7])U = 0, onfi, 

dV (6-2) 
— — = g on oil, 

ov 

Similarly, the Neumann function associated to the problem in the reference medium will be 
denoted by N w, (°). We denote by N w the Neumann function associated to the background 
medium, that is, 

(£ Ao , Mo +Po^ 2 [l + 7(x)])N"(x, y ) = -<5 y (x)I 2 , xett, x^y, 

<9N" ( 6 - 3 ) 
— — x,y =0 xeffl. 

ov 

We assume that 7 has small amplitude so that the Born approximation is valid. In particular, 
we have 

N"(x, y )c.N^(°)(x, y )+p oW 2 f N^ >(x,z)7(z)N^°)(z,y)dy. (6.4) 



As a consequence, we also have that U ~ — where 

U« (x) = -p co 2 f N W 'W (x, z) 7 (z)U( 0) (z)dz. (6.5) 



Let tr 7 denotes the typical size of 7, the remainders in the above approximations are of order 

o(cr 7 ). 

6.2 Statistics of the speckle field in the case of a density contrast only 

We assume that the inclusion has density contrast only. The backpropagation step constructs 
w as follows: 



w(x)=f r^(x,z)(i/-^' (0) )[U(")-u moas ](z)d t 7(z), xeO. (6.6) 
Jan z 

We emphasize that the backpropagation step uses the reference fundamental solutions, and the 
differential measurement is with respect to the reference solution. These are necessary steps 
because of the fluctuation in the background medium or equivalently, because of the fact that 
the background solution is unknown. 

Writing the difference between U'- -' and u mcas as the sum 

of u(°) - U and U - u mcas . These 
two differences arc estimated by U 1 - 1 ' in (|6.5[) and by (|2.22[) . respectively. Using Lemma [2.11 
we find that 



w(x) = - pquj 2 
-CS d 



[ rft(x-z)/ rS(z-y)U(°)(y) 7 (y)dyda(z) 
Jan Jn 

f r^(x - z)r|(z - z a )OT(z a )rfa(z) + OK^) + (a 7 ), xeft, 
Jan 



(6.7) 
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where C = uj 2 (po — pi)\B\. The second term is the leading contribution of U — u meas given by 
approximating the unknown Neumann function and the background solution by those associated 
to the reference medium. The leading error in this approximation is of order 0(a 1 8 d ) and can 
be written explicitly as 

C Po uj 2 5 d f T-(x,z) f rY(z,y)N^W(y,z Q )UW( Za ) 7 ( y )d y d (T (z) 
Jan Jn 

-Cp»LoH d f r£(x,z)rY(z,z a ) f N^W(z a ,y)UW( y ) 7 ( y )d yc ? ( T(z), 
Jan Jn 

and is neglected in the sequel. 

For the Helmholtz decomposition w Q , a £ {P, S}, the first fundamental solution Tq (x — z) 
in the expression (|6.7[) should be changed to Tq q (x — z). We observe that the second term in 
(|6.7[) is exactly (|3.10p . Therefore, we call this term w true and refer to the other term in the 
expression as w no ; sc . Using the Hclmholtz-Kirchhoff identity, we obtain 

w^ oisc (x) ~ f 7 (y)Q?m{r£ a (x - y)}Ul(y)dy, x e fl. (6.8) 

We have decomposed the backpropagation w" into the "true" w" ruo which behaves like in 
reference medium and the error part w" oise . In the TD imaging functional using multiple plane 
waves with equi-distributed directions, the contribution of w" rue is exactly as the one analyzed 
in Proposition 14.11 The contribution of w" oiso is a speckle field. 

The covariancc function of this speckle field, or equivalently that of the TD image corrupted 
by noise, is 

Cov(X W f[{U»}]( Z ),Z wp [{UJ}](z')) = C£-5 E^e{uf^(z).w^ nois Jz)}3?e{U i (0) ' Q (z).w« noisc (z)}], 

" 3,1=1 

for z, z' e fl, where C a is defined to be c a u> 2 \B'\(p' 1 / po — 1). Here U' )'" are the reference 
incoming plane waves ()3.14j) . 

Using the expression (|6.8|) . we have 




r 1 

= -b a / 7(y)- E eJK ° (z " y) ' eeje ^ ® e l ■ 5 ™{ r o> - yWy- 

where b a = (p Q uj)/c a . Finally, using ()3.18[) and (|3.19j) for a = P and S respectively, we obtain 
that 

1 " f 

-^Uf' tt (z).w« n0isc (z)^l / 7 (y)|3m{r^ Q (z-y)}| 2 rfy. (6.9) 
n ~pl Jn 

Here b' a = 4b a po(-^-) d ~ 2 (f^-) 2 ■ Note that the sum above is a real quantity. 
The covariance function of the TD image simplifies to 

Clb' a 2 I I C 7 (y,y')|3m{r^(z-y)}| 2 |3m{r^(z'- y ')}| 2 dydy', (6.10) 
JnJn 

where C 7 (y, y') = E[/y(y)7(y')] is the two-point correlation function of the fluctuations in the 
density parameter. 
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Remark 6.1. The expression in (|6.9|) shows that the speckle field in the image is essentially the 
medium noise smoothed by an integral kernel of the form I^suiTq a \ 2 . Similarly, f|6 . 10[) shows 
that the correlation structure of the speckle field is essentially that of the medium noise smoothed 
by the same kernel. Because the typical width of this kernel is about half the wavelength, the 
correlation length of the speckle field is roughly the maximum between the correlation length of 
medium noise and the wavelength. 



6.3 Statistics of the speckle field in the case of an elasticity contrast 

The case of elasticity contrast can be considered similarly. The covariance function of the TD 
image is 

n 

c 2 - E[5?e{VUf- a (z) : M'Vw« nolsc (z)}Se{VU} 0) ' a (z') : M'Vw£ noiso (z')}]. 

U 3,1=1 

Using the expression of w" oise , we have 

-£vU?(z) : M'Vw« noisc (z) = -b a / 7 (y)-E^« elK ° (Z " y) ' eej 
j=i JU 3=1 

eg } ® e a 6] ® el : [M'3m{V z r£ Q (z - y)}]dy. 

From (|3~T5|) and ([335]) . we see that 
1 " 

-^ a e^ e °ie gj ® e% ® e% = -4 W (-)^ 2 (^) 2 3m{V^ a (x)}. (6.11) 

I L f\ It, Qi 

3 = 1 

Using this formula, we get 

n 

-E VU iW : M 'Vw£ noisc (z) = b' a / 7 (y)Q 2 [M'](z - y)dy, (6.12) 

where Q 2 [M'](x) is a non-negative function defined as 

Q|[M](x) = 3m{Vr£ P (x)} : [M3m{Vr£ P (x)}] 

= a|3m{Vr^ P (x)}| 2 + &|3m{V • r^ P (x)}| 2 . 

The last equality follows from the expression (|2.19[) of M and the fact that diiY^ P )jk = 
dj(To p)ik- This symmetry is not satisfied for Tq p for which we have 

Q 2 [M](x) = 3m{Vr^ s (x)} : [MSm{Vr^ s (x)}] 

= ^|3m{Vr- 5 (x)}| 2 + |sSm{Vrgf iS (x)} : 3m{V^ s (x)} + 6|3m{V ■ T^ s (x)}| 2 . (6 ' 14) 

Here (V ■ Tq s (x))jki = dk(T^ s (x))ji- Note that Q 2 is non-negative and (|6.12l) is real valued. 
The covariance function of the TD image simplifies to 

[ [ C 1 (y,y')Ql[M'}(z-y)Ql[M\(z'-y')dydy', z,z'eO. (6.15) 
Jn Ja 

Remark 6.2. If we compare (|6.12|) with (|6.9j) . then one can see that they are of the same form 
except that the integral kernel is now Q 2 [M']. Therefore, Remark \6.1\ applies here as well. We 
remark also that the further reduction of the effect of measurement noise with rate 1/V2n does 
not appear in the medium noise case. In this sense, TD imaging is less stable with respect to 
medium noise. 
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6.4 Random elastic medium 



In this section we consider the case when the random fluctuation occurs in the elastic coeffi- 
cients. This is a more delicate case because it is well known that inhomogeneity in the Lame 
coefficients cause mode conversion. Nevertheless, we demonstrate below that as long as the 
random fluctuation is weak so that the Born approximation is valid, the imaging functional we 
proposed remains stable. 

To simplify the presentation, we assume that random fluctuation occurs only in the shear 
modulus fi while the density po and the first Lame coefficient Ao are homogeneous. The equation 
for time-harmonic elastic wave is then 

p a u 2 u + A V(V • u) + V • [/i(x)(Vu + (Vu) T )] = 0, (6.16) 

with the same boundary condition as before. The inhomogeneous shear modulus is given by 

p(x) = fi Q +7(x,w), (6.17) 

where 7(x) is a random process modeling the fluctuation. 

Born approximation. The equation for elastic wave above can be written as 

Cx ,^u + Po lu 2 u = -V • [7(x)(Vu + (Vu) T )]. 

Assume that the random fluctuation 7 is small enough so that the Born approximation is valid. 
We then have u ~ Uo — Ui where Uo solves the equation in the background medium and Ui , the 
first scattering, solves the above equation with u on the right hand side replaced by Uq. More 
precisely, we have 

m(x)= /N"(x,y)V-[7(y)(Vu + (Vu) T )]dy. (6.18) 
Jn 

Here, N" is the Neumann function in the background medium without random fluctuation. 

Post-processing step. As seen before, the post-processing (|3.6[) is a critical step in our 
method. As discussed in section 16. 2[ even when the medium is random we have to use the 
reference Green's function and the reference solution associated to the homogeneous medium 
in this post-processing step. Following the analysis in section we see that as in (|6.7[) the 
function w contains two main contributions: Firstly, back-propagating the difference between 
the measurement and the background solution in the random medium but without inclusion 
contributes to the detection of inclusion. Secondly, back-propagating the difference between 
the background solution and the reference solution in the homogeneous medium amounts to a 
speckle pattern in the image. 

The first contribution corresponds to the case with exact data and is discussed in section 4. 
We focus on the second contribution which accounts for the statistical stability. This part of 
the post-processed function w has the expression 

w? oiso (z)=H Q [^(i/-/C^)u 1 ]-H Q [/ rff(z,y)[(±r-£8) / N-(.,x)v(x)dx]( y )da( y )] 
L Jan L Jn 

= f r^(z,y) f Fjf(y,x)V • [ 7 (x)(Vu + (V^(x)]dxdtr(y). 
Jan Jn 

In the second equality above, v(x) is a short-hand notation for the divergence term in the line 
below. We refer to this term as the first scattering source. Using the Hclmholtz-Kirchhoff 
identity again, we obtain that 

w Q , noisc (z) = — / 3m{r^(z, x)}V • [ 7 (x)(Vu + (Vu) T )(x)]dx (6.19) 
c«w Jn 
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Remark 6.3. Compare the above expression with that in l|6.8p . The first scattering source in 
(|6.8[) is exactly the incident wave in the case with density fluctuation but is more complicated in 
the case with elastic fluctuation, see (|6.20[) below. This shows that the Born approximation in 
an inhomogeneous medium indeed captures weak mode conversion. Nevertheless, (|6.19[) shows 
that our method, due to the Helmholtz-Kirchhoff identity and our proposal of using Helmholtz 
decomposition, extracts only the modes that are desired by the imaging functional. As we will 
see, this is crucial to the statistical stability of the imaging functional. 

The speckle field. The specific expression of the imaging function depends on the type 
of the inclusion and the type of probing planewaves. We first consider the case of a density 
inclusion. For a pressure wave U p = e 4Kpx ee , the first scattering source is 

v(x) = 2w P V ■ (j(x)e lKpX - ee e e ® eg) = 2m P (V 7 • e e )U p (x) - 24 7 (x)U p (x). (6.20) 

The speckle field in the imaging functional with a set of pressure waves {U p } is 

.2 /Pi IMD'I 1 ' 



JwF,no ise [{U;}](z) = cp, 2 (^ - 1)|B'|-Ke£uf (z) • w p isc (z) 

3=1 

«(£i -l)|B'|»e / -24 7 (x)-y e ^ (z - x) - e9 , eej ®e e3 : 3m{r£ P (z, x)} 
Po Jn n ~[ 

n 

- 2- V lK p e iKp(z - x) -^ e 9j ® e 6j ® e 0j : [V 7 (x) ® 3m{^ P (z, x)}]dx. 



3=1 



Using the summation formulas (|3 . 1 8[) and (|6.11[) . we can rewrite the above quantity as 

Cf / 24 7 (x)|3m{r^ P (z,x)}| 2 +29m{V z r^(z,x)} : [V 7 (x) ® 3™{r£ P (z, x)}]dx. 
./fi 

Here, the constant is 



Cf = Vo (^- 2 (^)M^ - 1)1*1 = 4*"- a 4(Pi - Po)|B'| 



Assuming that 7 = near the boundary and using the divergence theorem, we can further 
simplify the expression of the speckle field to 

C[ [ 7 (x)[(24/ + A x )|3m{r^(z,x)}| 2 ]dx = C 1 p / [{2k 2 p I + A) 7 (x)]|3m{r^ P (z, x)}| 2 dx. 
Jn Jn 

This expression is again of the form of (|6.9[) and (|6.12[) except that the integral kernel is more 
complicated. Its correlation can be similarly calculated. Furthermore, Remarks 16.11 and 16.31 
apply here. The salient feature of the speckle field does not change. It is essentially the 
medium noise (or the derivative of the medium noise) smoothed by an integral kernel whose 
width is of the order of the wavelength. The correlation length of the speckle field is of the 
order of the maximum between that of the medium noise and the wavelength. 

The case when shear waves are used in (|5 . 1 1) can be similarly considered. For a shear wave 
U = e lK " sx ' eB e g L , the first scattering source becomes 

v(x) = iK S V ■ ( 7 (x)e- sx ' e « [eg ® ej + ej ® eg]) 

= -4 7 (x)U s (x) + ik s {V 1 ■ e e )U s (x) + m s (V 7 ■ ejy^^eg. 
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The speckle field in the imaging functional with a set of shear waves {U^} is 




- " E ^se iKs ^-^ e 9] ® ® e£ : [V 7 (x) ® 3m{r^(z, x)}] 
" i=i 

n 

_ -^i Ks e tes ( z - x > e «ie e x . <g> e£ <g> e*, : [V 7 (x) <g> Sm{r^(z, x)}]dx 
71 i=i 

Using the summation formulas (|3.19|l and (|6.11[) . we can rewrite the above quantity as 

Cf / 4 7 (x)|3m{r- s (z,x)}| 2 +^{V z r^ s (z,x)} : [V 7 (x) 8 3m{r^(z,x)}] 
+ Srn{(V z r^ s ) T (z,x)} : [V 7 (x) <g> 3m{r^(z, x)}]dx. 
Here, the constant is 

Cf = 4 Mo (-) d -M^ - = ^"^(P'x - Po)\B'\. 

K S PO Kg 

Moreover, the notation (VITq S ) T means the following: 

(vrg' iS (x))T w = ft(rs' i5 (x)) Jfc . 

Again, the speckle field is a smoothed version of the medium noise and its derivatives. The 
smoothing kernel can be read from the expression above. 

Elastic inclusion with shear waves. Now we consider the case of an elastic inclusion. 
The imaging function takes the form of (|5.15|> . When a set of pressure waves {U^} are used, 
the speckle field in the imaging function can be calculated as follows. 

1 " 

2wF,noi S c[{Uf}](z) = c P -VleJ2 VUf (z) : M'(B')Vw^ oi5C (z) 

= _$R e / - 7 (x)-^ym P e i ^ ( "- x ^ie 9] ®e * . <g> e e . : M'(B')3m{V z r£ P (z, x)} 
2 " 

- - y (iKp) 2 e iKp ( z -*'- ee i e 6j ® e 9j ® e 9j ® e 9] : [M'(B')Sm{V z r^ P (z, x)} ® V 7 (x)]dx. 

Using the summation formulas (|6.12[) and (|3.30j) . we can write the above quantity as 

C 2 P / 2Kp 7 (x)3m{Vr^p(z,x)} : [M'(B')3m{V z ^ P (z,x)}] 
Jn 

+ 23m{VVr^ P (z,x)} : [M'(S')3m{V z r^ P (z, x) ® V 7 (x)}]dx. 
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Here, the constant has the expression 

OJ Kp Kp Kp 

The tensor products above are defined as 

(M'(B')Sm{V z r% iP (z,x)}) jkl = m jkpq d p {Qm{T^p(z,x)}) q i, 
(M'(B')^{V z r^p(z,x) <g> V7(x)) iWa = (M'(B / )9m{V«rSf iP (z I x)}) iW a,7, 

(vvr^(z,x)) ifcl8 = a,a fc (r^ P (z,x)) ;s . 

Similarly, when a set of shear plane waves {U^} are used, the speckle field in the imaging 
function can be calculated as follows. 



1 " 

X W F,noi S e[{Uf}](z) = cs-SRe^VUf (z) : M'(S')Vwf oisc (z) 
n i=i 

i r k 2 " 

-Ke / - 7 (x)^^ lKs e lKs(z " x) - ee ,e ej ® e£ ® e£ : M'(B')MVz^ s (z, x)} 

n 

-^( lKs ) 2 e IKs(z - x) - efl ,e ej ®e ej . ® e£. ® e£. : [V 7 (x) ® M'(B')3m{V z r^ s (z, x)}] 



n 

- -y(iK S ) 2 e lKs{z -^e ej ®ei ®e e , ® e£ : [M'(B')3m{V z r^ P ( Z , x)} ® V 7 (x)]dx. 

j'=l 

Use formula (|6 . 1 2[) for the first integral and the differentiation of (|6.12[) for the second integral; 
use formula (|3.33[) for the third integral. We conclude that the speckle filed admits the following 
expression: 

Cl [ 4 7 (x)3m{Vr- s (z,x)} : ^(B')9m{V z r^(z,x)}] 

+ 3?m{VVr^ s (z,x)} : [V 7 (x) ® M'{B')Qm{V z r^ P (z, x)}] 
+ 3m{V 2 r^ s (z,x)} : [V 7 (x) ® M'(B')^m{V z T^ P (z, x)}]dx. 

Here, the constant has the expression 

C 2 s ^%o(-) d - 2 =4^- 2 ^. 

U) KS Kg 

Note also that the tensor WI^ s has coordinates {didj^T^ g)ki} and it is different from V 2 Tq S 
which is defined first in (|3 . 33|) . 

Remark 6.4. To summarize, we derived expressions for the speckle field in the imaging func- 
tional when the elastic medium has random fluctuations. We showed that the speckle field is 
essentially the medium noise (and its derivatives) smoothed by some integral kernels that de- 
pend only on the fundamental solution in the homogeneous medium. As a consequence, the 
covariance function of the speckle filed can be expressed as that of the medium noise (and those 
of its derivatives) smoothed by similar integral kernels. 
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7 Conclusion 



In this paper, wc performed an analysis of the topological derivative (TD) based clastic inclusion 
detection algorithm. We have seen that the standard TD based imaging functional may not at- 
tain its maximum at the location of the inclusion. Moreover, we have shown that its resolution 
does not reach the diffraction limit and identified the responsible terms, that are associated to 
the coupling of different wave-modes. In order to enhance resolution to its optimum, we can- 
celled out these coupling terms by means of a Hclmholtz decomposition and thereby designing 
a weighted imaging functional. We proved that the modified functional behaves like the square 
of the imaginary part of a pressure or a shear Green function, depending upon the choice of 
the incident wave, and then attains its maximum at the true location of the inclusion with a 
Rayleigh resolution limit, that is, of the order of half a wavelength. Finally, we have shown 
that the proposed imaging functionals are very stable with respect to measurement noise and 
moderately stable with respect to medium noise. 

In a forthcoming work, we intend to extend the results of the paper to the localization 
of the small infinitesimal clastic cracks and to the case of elastostatics. In this regard recent 
contributions [101 El [8] are expected to play a key role. 
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